##### ################################################### ######
#####                                                     ######
#####   Input: clean experiment data                      ######
#####   Output: power analysis                            ######
#####                                                     ######
##### ################################################### ######

setwd("/Users/lotte/Dropbox/PhD/style_experiment/replication")
rm(list=ls())

# Load libraries

library(data.table) # CRAN v1.14.2
library(plyr) # CRAN v1.8.6     
library(dplyr) # CRAN v1.0.9 
library(tidyverse) # CRAN v1.3.1
library(ggplot2) # CRAN v3.3.6 
library(paramtest) # CRAN v0.1.0

# Load clean experiment data

load("data/style_data.Rdata")
load("data/emotion_data.Rdata")
load("data/aggression_data.Rdata")
load("data/evidence_data.Rdata")

# Power analysis 

set.seed(190795)

lm_test_interaction <- function(simNum, N, b1, b2, b3, b0=0, x1m=0, x1sd=1,
                                x2m=0, x2sd=1) {
  
  ## Set x1 and x2 as binary
  x1 <- rbinom(N, 1, .5)
  x2 <- rbinom(N, 1, .5)

  y <- rnorm(N, b0 + b1*x1 + b2*x2 + b3*x1*x2, 1)
  model <- lm(y ~ x1 * x2)
  
  # Extract output from model (two main effects and interaction)
  est_x1 <- coef(summary(model))['x1', 'Estimate']
  p_x1 <- coef(summary(model))['x1', 'Pr(>|t|)']
  sig_x1 <- p_x1 < .05
  est_x2 <- coef(summary(model))['x2', 'Estimate']
  p_x2 <- coef(summary(model))['x2', 'Pr(>|t|)']
  sig_x2 <- p_x2 < .05
  est_int <- coef(summary(model))['x1:x2', 'Estimate']
  p_int <- coef(summary(model))['x1:x2', 'Pr(>|t|)']
  sig_int <- p_int < .05
  
  return(c(est_x1=est_x1, p_x1=p_x1, sig_x1=sig_x1, est_x2=est_x2, p_x2=p_x2,
           sig_x2=sig_x2, est_int=est_int, p_int=p_int, sig_int=sig_int))
}

# Assess power for b3 range(0.01,0.8)
power_lm_int <- grid_search(lm_test_interaction, params=list(b3 = seq(0.01, 0.8, 0.025)),
                            n.iter=5000, output='data.frame', b1=.2, b2=.2,  N=1650,
                            parallel='snow', ncpus=4)

results <- results(power_lm_int) %>%
  group_by(b3.test) %>%
  summarise(
    power_x1=mean(sig_x1),
    power_x2=mean(sig_x2),
    power_int=mean(sig_int)) 

pdf("analysis/plots/figure_S7_power_analysis.pdf", 6,6)
ggplot(results, aes(b3.test,power_int)) +
  geom_line() +
  theme_bw() + 
  geom_hline(yintercept = 0.8, lty="dashed") +
  xlab("Standardised effect size (conditional relationship between style treatment and MP gender)") +
  ylim(0,1) +  ylab("Power") +
  theme(axis.text = element_text(size = 9),
        axis.title = element_text(size = 9)) 
dev.off()
